## ctl+shift+n makes a new rscript

if (!require(pacman)) install.packages("pacman")
## Loading required package: pacman
pacman::p_load(tidyverse,
               sf,
               mapview,
               here)


# read/export vector data --------------------------------------------------

# read a shapefile (e.g., ESRI Shapefile format)
# `quiet = TRUE` just for cleaner output
sf_nc_county <- st_read(dsn = here("data/nc.shp"),
                         quiet = TRUE)

#save as shape file
st_write(sf_nc_county,
         dsn = here("data/sf_nc_county.shp"),
         append = FALSE)
## Deleting layer `sf_nc_county' using driver `ESRI Shapefile'
## Writing layer `sf_nc_county' to data source 
##   `/Users/mjparkersmi/Library/CloudStorage/OneDrive-UNCG/Documents/GitHub/cw_rgis/data/sf_nc_county.shp' using driver `ESRI Shapefile'
## Writing 100 features with 1 fields and geometry type Multi Polygon.
#save as geopackage
st_write(sf_nc_county,
         dsn = here("data/sf_nc_county.gpkg"),
         append = FALSE)
## Deleting layer `sf_nc_county' using driver `GPKG'
## Writing layer `sf_nc_county' to data source 
##   `/Users/mjparkersmi/Library/CloudStorage/OneDrive-UNCG/Documents/GitHub/cw_rgis/data/sf_nc_county.gpkg' using driver `GPKG'
## Writing 100 features with 1 fields and geometry type Multi Polygon.
#save as rds
saveRDS(sf_nc_county,
         file = here("data/sf_nc_county.rds"))

#read rds
sf_nc_county <- readRDS(file = here("data/sf_nc_county.rds"))


##point data 
sf_site <- readRDS(file = here("data/sf_finsync_nc.rds"))


mapview(sf_site,
        col.regions = "black", # point's fill color
        legend = FALSE) # disable legend
##take first 10 sites
sf_site10<- sf_site %>% slice(1:10)

mapview(sf_site10,
        col.regions= "black",
        legend = FALSE)
#line data
sf_str <- readRDS(file = here("data/sf_stream_gi.rds"))

mapview(sf_str,
        color = "steelblue", # line's color
        legend = FALSE) # disable legend
##take first 10 sites
sf_str10 <- sf_str %>% slice(1:10)

mapview(sf_str10,
        color = "darkblue",
        legend = FALSE)
#polygon
mapview(sf_nc_county,
        color="white",
        col.regions="orchid",
        legend= FALSE)
#highlight guilford county
sf_nc_gi <- sf_nc_county %>% 
  filter(county == "guilford")

mapview(sf_nc_gi,
        col.regions = "chartreuse4",
        legend = FALSE)
#map with ggplot
ggplot() +
  geom_sf(data = sf_nc_county) +
  geom_sf(data = sf_str) +
  geom_sf(data = sf_site) ##Not a good map :(

## a little better
ggplot() +
  geom_sf(data = sf_nc_gi) +
  geom_sf(data = sf_str)

# Exercises ---------------------------------------------------------------

#1.Read stream line data for Ashe county 
##Load the stream line data file sf_stream_as.rds located in the data folder. Use readRDS().
 sf_str_as <- readRDS(file = here("data/sf_stream_as.rds"))


#2. Check coordinate reference systems (CRS) (ref: Section 3.2.2)
##Print objects tocheck the CRS for both sf_str_as and sf_nc_county.
##Determine whether they share the same CRS.
print(sf_str_as)
## Simple feature collection with 6361 features and 1 field
## Geometry type: LINESTRING
## Dimension:     XY
## Bounding box:  xmin: -81.73891 ymin: 36.24445 xmax: -81.24562 ymax: 36.588
## Geodetic CRS:  WGS 84
## First 10 features:
##                          geometry       fid
## 1  LINESTRING (-81.32158 36.47... fid000001
## 2  LINESTRING (-81.31133 36.47... fid000002
## 3  LINESTRING (-81.38594 36.45... fid000003
## 4  LINESTRING (-81.39307 36.46... fid000004
## 5  LINESTRING (-81.3454 36.533... fid000005
## 6  LINESTRING (-81.3472 36.532... fid000006
## 7  LINESTRING (-81.3328 36.518... fid000007
## 8  LINESTRING (-81.34964 36.50... fid000008
## 9  LINESTRING (-81.35949 36.53... fid000009
## 10 LINESTRING (-81.35094 36.52... fid000010
print(sf_nc_county)
## Simple feature collection with 100 features and 1 field
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -84.32377 ymin: 33.88212 xmax: -75.45662 ymax: 36.58973
## Geodetic CRS:  WGS 84
## First 10 features:
##         county                       geometry
## 1         ashe MULTIPOLYGON (((-81.47258 3...
## 2    alleghany MULTIPOLYGON (((-81.23971 3...
## 3        surry MULTIPOLYGON (((-80.45614 3...
## 4    currituck MULTIPOLYGON (((-76.00863 3...
## 5  northampton MULTIPOLYGON (((-77.21736 3...
## 6     hertford MULTIPOLYGON (((-76.74474 3...
## 7       camden MULTIPOLYGON (((-76.00863 3...
## 8        gates MULTIPOLYGON (((-76.56218 3...
## 9       warren MULTIPOLYGON (((-78.30849 3...
## 10      stokes MULTIPOLYGON (((-80.02545 3...
##CRS for sf_str_as: WGS 84
##CRS for sf_nc_county: WGS 84

###geodectic CRS for both is WGS 84

#3. Export to GeoPackage format (ref: Section 3.2.2)
##Save sf_str_as as a GeoPackage file named sf_stream_as.gpkg in the data folder. Use sf::st_write().

###SKIP not working?
# st_write(sf_str_as, 
#          dsn = "data/sf_stream_as.gpkg",
#          append = FALSE)


#4. Map streams and county boundaries (ref: Section 3.2.6)
##Create a map displaying both: North Carolina county boundaries from sf_nc_county, 
##and Ashe County stream lines from sf_str_as. Use ggplot2::ggplot() with multiple geom_sf() layers.

ggplot() +
  geom_sf(data = sf_nc_county) +
  geom_sf(data = sf_str_as) ##not a good map again :(

#5. Subset county layer to Ashe county and remap (ref: Section 3.2.5 and Section 3.2.6)
##Subset the sf_nc_county object to include only Ashe County. Use dplyr::filter() for subsetting. Assign the result to sf_nc_as.
##Then, recreate the map showing only sf_nc_as and sf_str_as
sf_nc_as <- sf_nc_county %>%
  filter(county == "ashe")

ggplot() +
  geom_sf(data = sf_nc_as) +
  geom_sf(data = sf_str_as, color = "dodgerblue")